Phylogeography and cohesion species delimitation of California endemic trapdoor spiders within the Aptostichus icenoglei sibling species complex (Araneae: Mygalomorphae: Euctenizidae)

Abstract Species delimitation is an imperative first step toward understanding Earth's biodiversity, yet what constitutes a species and the relative importance of the various processes by which new species arise continue to be debatable. Species delimitation in spiders has traditionally used morphological characters; however, certain mygalomorph spiders exhibit morphological homogeneity despite long periods of population‐level isolation, absence of gene flow, and consequent high degrees of molecular divergence. Studies have shown strong geographic structuring and significant genetic divergence among several species complexes within the trapdoor spider genus Aptostichus, most of which are restricted to the California Floristic Province (CAFP) biodiversity hotspot. Specifically, the Aptostichus icenoglei complex, which comprises the three sibling species, A. barackobamai, A. isabella, and A. icenoglei, exhibits evidence of cryptic mitochondrial DNA diversity throughout their ranges in Northern, Central, and Southern California. Our study aimed to explicitly test species hypotheses within this assemblage by implementing a cohesion species‐based approach. We used genomic‐scale data (ultraconserved elements, UCEs) to first evaluate genetic exchangeability and then assessed ecological interchangeability of genetic lineages. Biogeographical analysis was used to assess the likelihood of dispersal versus vicariance events that may have influenced speciation pattern and process across the CAFP's complex geologic and topographic landscape. Considering the lack of congruence across data types and analyses, we take a more conservative approach by retaining species boundaries within A. icenoglei.


| INTRODUC TI ON
The conceptual definition of what constitutes a species along with the relative importance of the varied processes by which new species arise continue to be much-debated topics of discussion (de Queiroz, 2007;Hey, 2001;Wells et al., 2021). Species concepts typically emphasize disparate intrinsic biological properties (e.g., morphological differences, niche divergence, and genetic divergence) that can be differentially important with respect to species recognition and/or speciation process. Contingent factors, for example taxon characteristics and life history traits and point/stage in the speciation process, may render various concepts incompatible and/or delimit species in different ways; that is, one concept may recognize multiple distinct species whereas another may lump them together (de Queiroz, 2007).
Assessing species limits is particularly difficult in taxa with limited dispersal capabilities when reduced gene flow leads to high levels of population structuring. Taxa with high levels of genetic divergence and no gene flow can sometimes lead to speciation in the absence of notable morphological differentiation, essentially obscuring species boundaries. Specifically, non-vagile taxa are closely tied to the landscape, such that as geological, topographical, and climatic changes occur over time, populations become geographically isolated with severely limited opportunity for gene flow (Bond et al., 2001;Derkarabetian et al., 2021;Starrett & Hedin, 2007;Weisrock & Larson, 2006). As these populations remain spatially isolated over relatively long periods of time and accumulate random mutations, genetic divergence builds through genetic drift and/or natural selection for adaptive alleles in population(s) that inhabit newly available niche space. When spatial isolation is coupled with occupation of new niche space, one would expect each population to not only exhibit genetic divergence but also morphological, behavioral, and/ or physiological differences (Freudenstein et al., 2016). However, if genetically diverged populations remain stationary in niche space (i.e., niche conservatism; Wiens & Graham, 2005), then it would be plausible for morphological stasis to occur in the absence of differing selective pressures, resulting in genetic lineages that are morphologically indistinguishable (Bond et al., 2001;Cerca et al., 2021;Mas-Peinado et al., 2018). In that case, it is likely that species diversity will be underestimated because traditional approaches that primarily apply morphological distinctiveness are more commonly used in species delimitation (Bond et al., 2021). Thus, implementing a species concept focusing on one biological property/data type could potentially misrepresent the actual number of species present if that property was not important in the speciation process (Abbott et al., 2013;de Queiroz, 2007).
The species concept applied in a given system has implications for downstream delimitation decisions and thus must be explicitly stated in any species delimitation study. Nevertheless, in many taxonomic studies (e.g., in spider taxonomy), an explicit species concept is seldom stated (Bond et al., 2021). A species concept that focuses strictly on morphological differentiation has the potential to overlook cryptic species that may otherwise be genetically diverged to the point that genomic incompatibilities preclude gene flow (Barroso et al., 2010;Battey & Klicka, 2017;Holland et al., 2004;Weisrock & Larson, 2006). Alternatively, molecular approaches to species delimitation have been shown to overestimate species diversity when local population structuring is viewed as "species divergence".
Specifically, single-locus approaches such as DNA barcoding along with GMYC, as well as multiple-locus approaches (e.g., multispecies coalescent methods) that assume panmixia are prone to identifying population structuring as opposed to speciation events (Hamilton et al., 2014;Hedin et al., 2015;Sukumaran & Knowles, 2017). In such systems, localized divergence of neutral alleles may be inconsequential when populations come into secondary contact, so any species delimitation approach that relies primarily on genetic differentiation has the potential to overestimate species diversity when applying these methods. Consequently, a species concept that incorporates multiple biological properties as an integrative species delimitation approach that weighs evidence from multiple independent sources is likely to more accurately identify true evolutionary species diversity.
The Cohesion Species Concept (CSC; Templeton, 1989) has arguably already solved the problems of too little versus too much gene flow and provides the hypothetical and conceptual foundation for framing integrative species delimitation. The CSC posits that a cohesion species must constitute an independently evolving evolutionary lineage and must be genetically exchangeable and/or ecologically interchangeable (Templeton, 1989). Specifically, a primary tenet of a cohesion species is that it comprises populations that exchange genes and occupy similar niche space. This concept can be applied to essentially all taxa, integrates multiple biological properties that are potentially important in the speciation process, and provides a methodological framework in which species hypotheses can be tested (Barraclough, 2019;Bond & Stockman, 2008;Templeton, 1989;Wells et al., 2021). Thus, it is particularly useful when evaluating species boundaries in morphologically homogenous taxa prone to cryptic diversity in conjunction with a high amount of population structuring at small spatial scales (Bond & Stockman, 2008;Hendrixson et al., 2013Hendrixson et al., , 2015Newton et al., 2020).
In this paper, we will apply the CSC to a species delimitation problem in a previously characterized lineage of trapdoor spiders in the genus Aptostichus Simon (Araneae: Mygalomorphae: Euctenizidae), specifically species in the Aptostichus icenoglei sibling species complex. Mygalomorph spiders are notorious for being morphologically static relative to the other, more diverse spider groups placed in the sister infraorder Araneomorphae (Opatova et al., 2019); they have relatively long lifespans and limited dispersal capabilities which makes their populations more vulnerable to genetic structuring at small spatial scales (Bond & Stockman, 2008;Cooper et al., 2011;Hamilton et al., 2014;Harvey et al., 2018;Starrett & Hedin, 2007), thus underscoring the interplay of genetic versus ecological interchangeability when evaluating divergence at the species/population interface in these highly structured taxa. The questions we pose are related first to genetic exchangeability-do these populations constitute distinct genetic lineages, and if so, are they ecologically interchangeable, or not? If these genetically distinct lineages are ecologically interchangeable, then the unavoidably subjective question arises of how heavily one weights genetic divergence versus ecological/ adaptive divergence, or the lack thereof (discussed below).
Aptostichus barackobamai and A. icenoglei are relatively widespread and exhibit evidence of cryptic diversity (i.e., morphologically similar yet found in a variety of habitats across a sizable geographic range) found in other mygalomorph groups (Hamilton et al., 2014;Hendrixson et al., 2013;Hendrixson & Bond, 2005;Starrett et al., 2018;Starrett & Hedin, 2007) as well as other Aptostichus species (Bond et al., 2001;Bond & Stockman, 2008). Aptostichus isabella, on the contrary, is only known from one specimen collected near Lake Isabella in the southern Sierran foothills. Aptostichus icenoglei is distributed throughout the Transverse Ranges from the Los Angeles Basin to the Santa Ana/San Jacinto Mountains as well as the mountains and hills surrounding San Diego (Bond, 2012). The primary habitat types for A. icenoglei include coastal chaparral forest and coastal range open woodland shrub and coniferous forest (Bond, 2012). Aptostichus barackobamai is found in primarily mixed redwood and coniferous forests in the northern Coastal Ranges as well as along the northern rim of the Central Valley, with one population in the Sutter Buttes (Bond, 2012). Altogether, these likely represent a diversity of habitat types spread across a number of different California ecoregions. Mitochondrial data from Bond (2012) indicate considerable population genetic structuring, especially in A. icenoglei, which is likely influenced by the typical mygalomorph life history traits discussed above. This, in conjunction with notable molecular divergence as well as a diversity of habitats, suggests that A. barackobamai and A. icenoglei populations, respectively, have been isolated from gene flow for an extended period of time, which would increase speciation potential (i.e., both likely comprise more than one species; Barraclough, 2019).
The primary objective of this study was to use multiple lines of evidence, specifically morphological, ecological, and genomic-scale data (i.e., ultraconserved elements, UCEs; Faircloth et al., 2012) and to evaluate phylogenetic relationships, species boundaries, and historical biogeography within the A. icenoglei complex. We explicitly tested species hypotheses within this assemblage by implementing a CSC-based approach. We first evaluated genetic exchangeability using clustering analyses to assess the potential for gene flow and then assessed ecological interchangeability of genetic lineages with a niche-based distribution modeling approach. Additionally, biogeographic analysis was used to investigate the likelihood of dispersal versus vicariant events that may have influenced speciation pattern and process across the CAFP's complex topographic and geologic landscape.

| Taxon sampling
We sampled 62 individuals overall for the three species within the complex using both specimens from Bond (2012) and new records ( Figure 1; see Table S1 for locality information). Aptostichus barackobamai was collected across its geographic range in northern California for a total of 21 samples, and A. icenoglei was collected throughout its range in southern California for a total of 40 samples. Only one specimen of A. isabella was included in this study due to collecting constraints (i.e., only one individual of this species has ever been collected and a burrow has not yet been found containing this species; Bond, 2012).

| Sequence capture
Data for ultraconserved elements were produced following the Sequence processing and analyses were performed on the Farm Community HPC at the University of California, Davis. Reads were filtered and trimmed using Illumiprocessor (Faircloth, 2013) and Trimmomatic (Bolger et al., 2014) in the Phyluce 1.7.1 pipeline (Faircloth, 2015). De novo assemblies with the cleaned paired-end and single-end reads were performed using SPAdes v. 3.14.1 with the isolate option (Prjibelski et al., 2020). Scaffolds were matched with 65% identity and 65% coverage to the modified probe list from Maddison et al. (2020), which is a blend of the Arachnid (Faircloth, 2017;Starrett et al., 2017) and Spider (Kulkarni et al., 2020) probe sets. MAFFT (Katoh & Standley, 2013) was used to align individual locus datasets, and alignments with locus occupancy (i.e., completeness) minimums of 50%, 75%, and 90% were obtained. Alignment masking was performed with TrimAl v.1.2 (Capella-Gutierrez et al., 2009) using default settings.
SNP datasets were generated for A. icenoglei only with Phyluce from the 50%, 75%, and 90% minimum occupancy loci. Reads were mapped against corresponding scaffolds with BWA (Li & Durbin, 2009), implemented in Phyluce, and phased alignments were generated for each minimum locus completeness set. Phased alignments were screened for SNPs, with five sets of single random SNP per locus generated to test for SNP set sensitivity.

| Phylogenetic & biogeographic analyses
Phylogenies were estimated for three different data sets (50, 75, and 90 percent locus completeness; Figure 2 and Figures S1 and S2) with a maximum likelihood inference using IQ-TREE v2.1.2 (Minh et al., 2020). Model selection was performed by ModelFinder (Kalyaanamoorthy et al., 2017), which is implemented in IQ-TREE, and support values were inferred from 1000 replicates of ultrafast bootstrapping (Hoang et al., 2018). Our phylogenies were visualized in FigTree v1.4.1 with midpoint rooting (midpoint rooting produces a result consistent in other analyses; Bond, 2012) and compared to assess congruence among clades. We also conducted two coalescentbased analyses for the 75p and 90p data sets. Gene trees for each locus were constructed using RAxML v8.0.12 for each data set and used to generate a coalescent-based tree with ASTRAL-III (Zhang et al., 2018). Multispecies coalescent (MSC) bootstrapping was run with ASTRAL v.5.7.4 and 100 pseudoreplicates (Simmons et al., 2019). For downstream analyses, we employed the ML phylogeny based on the largest amount of taxon coverage and with the most robust support values (i.e., the phylogeny with 90 percent minimum locus completeness).
Biogeographic analyses were generated using Reconstruct Ancestral State in Phylogenies (RASP; Yu et al., 2015) with dispersal constraints (i.e., dispersal multipliers set to 0.01 for adjacent areas and 0.0001 for non-adjacent areas) to account for their limited dispersal capacity and using our 90p consensus tree from IQ-TREE.
Model testing was conducted using the R package BioGeoBEARS

F I G U R E 1 Geographic distributions of Aptostichus icenoglei
sibling species complex lineages. Inset legend denotes color scheme for each of the lineages recovered in Figure 2.

| Cohesion species delimitation
To assess species boundaries within A. barackobamai and A. icenoglei, we employed the methodological framework for delimiting cohesion species from Bond and Stockman (2008) that evaluates two cohesion mechanisms: genetic exchangeability and ecological interchangeability. We used our 90p topology as the baseline evolutionary framework for establishing the "basal starting point" to  within A. barackobamai, we designated all the individuals as part of one evolving lineage that was not tested further for genetic and ecological exchangeability. For A. icenoglei, we also used our topology from the MSC tree resampling ( Figure S6) as additional guidance for establishing lineage designations, which resulted in 3 lineages: North, Central, and South (see Figure 2). We evaluated the distributions of these lineages as well as performed morphological and genetic clustering analyses to assess the potential for gene flow.
Genetic exchangeability was rejected if any allopatric lineage forms an apparently separate clustering pattern from other lineages, or if any parapatric lineage has a separate clustering pattern and an obvious barrier to gene flow.
For morphological data, we quantified 25 continuous character measurements for 30 males (10 males from each lineage; Table S2).
All measurements were recorded in millimeters and were quantified with a Leica M165C stereomicroscope using the Leica Application Suite software and a digital camera. Measurements were trans- SNP datasets were converted to one-hot encoding, which converts categorical data into numerical data as needed for certain machine  Fick & Hijmans, 2017). Climate layers were cropped to encompass the geographic area of interest and converted to a raster stack using R packages raster (Hijmans, 2015) and rgdal (Bivand et al., 2019). Highly correlated variables with a Pearson correlation coefficient > .80, estimated using the R package ENMTools (Warren et al., 2021), were removed. The remaining bioclimatic variables (see Table S3) were used in conjunction with occurrence records from the current study as well as records from Bond (2012) (Phillips & Dudík, 2008), which is a machine learning program that uses a maximum entropy algorithm.
Multiple points within a 30 arc-second grid cell were removed (i.e., only retaining one record per grid cell) by ENMeval during the modeling step to reduce potential for spatial autocorrelation. To limit the likelihood of overfitting while also accounting for goodness of fit, multiple feature classes and regularization multipliers were chosen to generate a total of 30 models (see Tables S4-S7 for model parameters and stats). Model selection was based on AICc, with the best model having a delta AICc of zero and was subsequently used in downstream analyses ( Figure 5).
Statistical comparisons of SDMs for each sister lineage comparison were conducted with niche overlap, niche equivalency, and niche similarity tests in ENMTools (Warren et al., 2008(Warren et al., , 2010. We used the Schoener's D statistic (Schoener, 1968) to calculate the niche overlap for each lineage comparison, which ranges from 0 (no overlap) to 1 (complete overlap). We carried out two tests, niche equivalency and niche similarity, to evaluate the significance of niche overlap with a randomization procedure (Warren et al., 2008(Warren et al., , 2010. The niche equivalency test, a one-tailed test, assesses whether the two niches being compared are identical or not. If the observed niche overlap value is significantly lower than the null distribution of randomized D values, then the niches are not identical (i.e., not equivalent; Figure S7). Considering the limitation of relying only on occurrence records for the niche equivalency test (Warren et al., 2008), we also employed the niche similarity test, a two-tailed test, to assess whether niche overlap between lineages relative to the niche spaces available to those lineages is more similar or different than expected by chance (niche conservatism or (2) minimum area polygons based on SDM raster grid cells with a habitat suitability score threshold greater than 0.5 (i.e., a polygon was generated around every grid cell with a habitat suitability score greater than 0.5); and (3) minimum area polygons based on SDM raster grid cells with a habitat suitability score threshold greater than 0.75 (see Figures S10-S12 for reference).

| UCE stats
The UCE data are summarized in Table 1

| Phylogenetic and biogeographic analyses
All estimated phylogenies fully supported (i.e., 100 for IQ-TREE or one for ASTRAL analyses) species level divergence among the three previously delineated morphological species within this sibling complex (see Figure 2 and Figures S1-S6). Also, all A. icenoglei lineage   Figure 1). Three clustering analyses, one with morphological data and two with molecular data, were also used to inform the possibility of gene flow. PCA analysis of the quantitative morphological measurements reveal no distinct clustering for any of the lineages (Figure 4). Similarly, the CLADES analysis with the Metano_CLADES training model identified one species. In contrast, the VAE analysis for 50p indicates three very distinct clusters corresponding to each lineage for both the mean and standard deviation ( Figure S13); however, although VAE analyses for 75p and 90p

| Cohesion species delimitation
show separation between the lineages for the mean, there is a small amount of overlap for the standard deviation between Central and South lineages (Figure 4). South occurrence points to Central background was significantly more similar than expected by chance (i.e., niche conservatism; Figure S8). Central occurrences compared with the South background, determined by minimum bounding polygons based on raster grid cells with either habitat suitability scores >0.5 or >0.75, and vice versa indicated niche conservatism (i.e., more similar than expected compared with the null distribution; Figure S8). When comparing the Central+South occurrence records to the minimum bounding polygon connecting occurrence points defining the background region of North, the results show no significant difference; yet, the reciprocal comparison is significantly more similar than expected ( Figure S9). All comparisons of North versus Central+South and vice versa suggest niche conservatism when background regions are defined by minimum bounding polygons based on raster grid cells with either habitat suitability scores >0.5 or >0.75 ( Figure S9).

| DISCUSS ION
species-based delimitation approach from Bond and Stockman (2008) to expand on the evaluation of species boundaries within the complex.
Once evolutionary lineages were delineated, based on the topology and high support values (i.e., >0.95) from both the 90p IQ-TREE and the 90p MSC bootstrapping (Figure 2 and Figure

| Speciation and phylogeography
Spiders in the A. icenoglei complex, like most mygalomorphs, have limited dispersal capabilities and relatively long generation times (Bond, 2012;Harvey et al., 2018;Hedin et al., 2013;Hendrixson et al., 2013), which contributes to their tendency to have population structure at relatively small spatial scales. The molecular data show genetic divergence across the A. icenoglei complex and within A. icenoglei populations, thus populations have likely been isolated from gene flow for a long period of time, indicating the increased potential for speciation (Barraclough, 2019 (Hendrixson et al., 2013;McCormack et al., 2009;Newton et al., 2020;Starrett et al., 2018;Warren et al., 2008), yet very few are explicit about the background region they chose to incorporate into the analysis (McCormack et al., 2009;Newton et al., 2020;Starrett et al., 2018).
In addition, as far as we are aware, no other study other than Newton Ranges approximately 5 mya (Norris & Webb, 1990), which likely cut off the potential for gene flow, and has been hypothesized for other CAFP taxa (Alexander & Burns, 2006;Calsbeek et al., 2003;Feldman & Spicer, 2006;Reilly et al., 2015;Rissler et al., 2006). Second, the split between A. isabella and A. icenoglei possibly occurred due to vicariance. This split could be attributed to the Tehachapi Mountains acting as a barrier to dispersal, which has also been inferred for other taxa (Calsbeek et al., 2003;Chatzimanolis & Caterino, 2007;Rissler et al., 2006). Third, vicariance was inferred for the split between the North lineage and Central+South lineages, which could be associated with periodic inundations of the LA Basin (Jacobs et al., 2004) that might have resulted in habitat fragmentation, also hypothesized for the mahogany Jerusalem cricket (Stenopelmatus "mahogani"; Vandergast et al., 2006) and stream-dwelling frogs (Pseudacris cadaverina; Phillipsen & Metcalf, 2009).

| Species limits within Aptostichus icenoglei and taxonomic implications
Although our integrative approach considered multiple independent lines of evidence, our conflicting results circle back to the unavoidably subjective question of how much weight should be given TA B L E 2 Summary of Aptostichus icenoglei cohesion species delimitation assessment.

Lineage comparison
Genetic exchangeability

CLADES (molecules) Conclusion
Central to South Parapatric, no obvious barrier to genetic divergence versus morphological/ecological divergence (or lack thereof) when delimiting species with extreme population structuring. Should we elevate genetically diverged lineages to species status despite the lack of observed morphological/ecological differences? One could argue that identifying and describing cryptic diversity can be important not only for more accurate biodiversity measures but also conservation management plans (i.e., evolutionary significant units; Ryder, 1986). For example, Fennessy et al. Alternatively, one could argue that there is no practical value of recognizing genetically diverged lineages as separate species considering the lack of any visible diagnostic character/difference in ecological role (Freudenstein et al., 2016). Specifically, Freudenstein et al. (2016) argued that possessing both a unique ecological role and phenotypic differences are imperative when recognizing distinct species units. However, even this argument is rife with subjectivity; for example, how much phenotypic difference is enough to distinguish lineages as separate species? Also, it has been established that the speciation process is a continuum in which certain biological properties can be affected at different points along that continuum (Abbott et al., 2013;de Queiroz, 2007). Thus, it is feasible for geographically separated populations to accumulate enough Thus, the most conservative taxonomic approach would be to require rejection of both genetic and ecological interchangeability for identifying separate cohesion species.
Studies spanning different animal taxa that have utilized CSCbased delimitation approaches have highlighted the importance of evidence for adaptive divergence when delimiting species (Bond & Stockman, 2008;Leaché et al., 2009;Newton et al., 2020;Rengifo-Correa et al., 2021). For mygalomorphs, Bond and Stockman (2008), the study upon which our CSC framework is based, delimited Our analytical results separately inferred one to three species within A. icenoglei depending on the dataset and analysis used, but the final decision, arguably subjective, comes down to emphasizing mygalomorph life history characteristics and acknowledging limitations for each data type/analysis (discussed further below). The three species hypothesis was dismissed due to: (1) the less conservative niche equivalency test (Warren et al., 2008), (2) the possibility that the 50 percent locus completeness SNP dataset overly inflated cluster separation between Central + South, and (3) no obvious barrier to gene flow between Central and South lineages. The two species hypothesis is not substantiated based on morphological and ecological similarity between lineages, yet it is supported by rejecting genetic exchangeability as inferred by the VAE cluster separation with higher/more conservative locus completeness percentage datasets (75p and 90p) and a probable hard barrier to gene flow, the LA Basin, between North and Central+South. Although the LA Basin is likely impeding gene flow due to urban development and habitat fragmentation, the small likelihood of a potential corridor of habitat suitable for dispersal along the northern Basin rim/southeastern San Bernardino mountains cannot be completely dismissed (e.g., figure 1 in Vandergast et al., 2006). The one species hypothesis is supported by morphological and ecological data as well as an implementation of a supervised machine learning analysis on the 90p SNP dataset.
The flowchart in Bond and Stockman (2008)

| Limitations of analyses and future prospects
We believe that the supervised machine learning analysis has limitations due to potential shortcomings with the training data set devised using unrelated taxa. Although we see great value in attempting to establish a training dataset integrating biologically/ ecologically relevant characteristics, it is difficult to assess how applicable this dataset can be to other dispersal-limited taxa, especially across different taxonomic orders and biogeographical regions (Derkarabetian et al., 2022 (Myers et al., 2000).
One could argue that the overall complexity of the CAFP might influence the speciation process of low dispersal taxa in a different manner from how topographic changes in the Pacific Northwest would to the point that the genetic signatures may manifest differently.
Specifically, as there are more topographical changes (both in number and intensity), the more chances there are for speciation through vicariance when compared to fewer/less drastic topography shifts (Badgley et al., 2017).
Our VAE analysis with the lower locus completeness dataset (50p) showed obvious separation between all three of the A. icenoglei lineages, whereas our higher locus completeness datasets (75p and 90p) retained only enough signal to maintain the North lineage as a separate cluster but not for Central or South lineages (Figure 4). VAE relies on the inherent structure present in the data (Derkarabetian et al., 2019), and previous studies have shown that VAE analyses have been heavily influenced by the filtering parameters for the SNP datasets (Martin et al., 2021;Newton et al., 2020). Specifically, if a lower threshold for locus completeness is allowed in a dataset the more likely it is to "over-split", whereas more stringent filtering (i.e., a high threshold for locus completeness) can remove potentially important signal and "under-split" the amount of diversity. Because our higher locus completeness datasets retained the same clustering patterns, we are confident that they accurately reflect genetic divergence, and that the 50p dataset separation pattern for Central and South is an artifact of the filtering choice. Thus, it is important to be mindful of the potential filtering strategies for these SNP datasets, and best practices if utilizing VAE as a species delimitation method would be to use multiple filtering strategies to identify possible data artifacts versus actual structure.
There are known caveats for using niche-based distribution modeling approaches as a proxy for evaluating ecological interchangeability. First, it has to be acknowledged that large-scale ecological data, which are based on a very small time frame of 30 years (i.e., 1970-2000), used for building the SDMs potentially lacks the resolution needed for detecting very small-scale habitat differences which may be important for detecting adaptive divergence (Massatti & Knowles, 2014;Newton et al., 2020;Starrett et al., 2018). The microhabitat preferences for these spiders, which includes shaded ravines, north-facing slopes, and specific soil types (Bond, 2012), found within the heterogeneous landscapes throughout the CAFP are potentially not identified in the SDMs by even the best resolution available. Thus, our niche similarity tests using these models likely do not detect the potential for microhabitat niche divergence and consequently suggest the need for studies that assess fine-scale data on variables like temperature, precipitation, burrow features (e.g., size and depth), and soil composition.
Second, as discussed above, background region choice can heavily impact the results of niche similarity tests, thus incorporating multiple regions with biologically relevant constraints may provide a more rigorous application. Third, considering that our proxy of ecological interchangeability was only based on the abiotic factors contributing to niche space (i.e., bioclimatic variables and occurrence records to build an SDM), one could argue that there were other potential avenues of ecological divergence that could have been included in this study for a more robust evaluation of ecological interchangeability. There are potential biotic factors (e.g., competition with other taxa, difference in prey items across microhabitats, or non-overlapping breeding periods) that could distinguish lineages from one another. For example, other studies delimiting mygalomorph species have included behavioral traits when applicable (e.g., non-overlapping breeding periods; Hendrixson et al., 2015;Hendrixson & Bond, 2005;Prentice, 1997). Unfortunately, the lack of available natural history data for many fossorial mygalomorphs (Bond, 2012;Hedin et al., 2013;Starrett et al., 2017) have limited use of this type of data in species delimitation decisions.
Given these limitations, there are many potential avenues in which researchers can begin to bridge these gaps in knowledge.
First, generating more datasets comprising low-dispersal taxa with varying divergence times and across other biogeographical regions that can be used to train models for supervised machine learning methods such as CLADES, will likely aid the robustness of this approach (Derkarabetian et al., 2022). Second, accumulating more natural history data for mygalomorphs will not only provide valuable general ecological information but may also be used as additional evidence in species delimitation. For example, pitfall trapping spiders in areas where occurrence records of each species/lineage of interest is well-known to collect specimens can be informative for both breeding period times and gut content analysis to identify prey items that are being ingested (i.e., can inform potential for ecological divergence). Third, the advent of assembled and annotated genomes for non-model taxa, specifically in Aptostichus, will likely pave the way toward utilizing these data not only for reconstructing evolutionary relationships but also identifying genes that contribute to potential adaptive divergence across the landscape (Johnson, 2018).
Overall, our study emphasized the efficacy of implementing a cohesion species-based delimitation approach across all taxa, but especially for assessing the potential of cryptic diversity. Using genome-wide UCEs in conjunction with morphological and ecological data to evaluate genetic and ecological exchangeability provided multiple independent lines of evidence that covered multiple biological properties potentially important in the speciation process.
Specifically, this integrative approach underscored how different data types or approaches alone could either over-or under-split diversity estimates, yet taking them all into consideration led to a more robust species delimitation hypothesis within the A. icenoglei complex. Typically, such studies of taxa with extreme population structuring favor recognizing cryptic species, whereas herein, we have shown that an integrative approach, considering multiple lines of evidence, has the capacity to retain (lump) populations as a single species. Moreover, we reinforce the capability of the Cohesion Species Concept in providing both the conceptual and experimental framework for conducting such tests. Finally, our biogeographic analysis reveals that vicariance likely played a dominant role in the speciation process across the entire complex, further highlighting the impact of the complex geological, climatic, and topographical changes across the CAFP on speciation process.

CO N FLI C T O F I NTER E S T S TATEM ENT
Collectively, the authors (LGN, JS, EEJ, JEB) have no conflict of interest to disclose.

DATA AVA I L A B I L I T Y S TAT E M E N T
All sequence data can be found in Sequence Read Archive with the BioProject ID PRJNA949729. All tree files, data matrices, and VAE one-hot encoded files can be found here: https://doi.org/10.25338/ B8MD2Z. All scripts used in this study can be found here: https:// github.com/lgnew ton/AptIc eClade_SpDelim.